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The collective dynamics of coupled two-dimensional chaotic maps on complex networks is 
known to exhibit a rich variety of emergent properties which crucially depend on the underlying 
network topology. We investigate the collective motion of Chirikov standard maps interacting 
with time delay through directed links of Gene Regulatory Network of bacterium Escherichia Coli. 
Departures from strongly chaotic behavior of the isolated maps are studied in relation to different 
coupling forms and strengths. At smaller coupling intensities the network induces stable and 
coherent emergent dynamics. The unstable behavior appearing with increase of coupling strength 
remains confined within a connected sub-network. For the appropriate coupling, network exhibits 
statistically robust self-organized dynamics in a weakly chaotic regime. 

The system of transcriptional regulations among 
genes in a cell is suitably modeled by Gene Reg- 
ulatory Network, with nodes representing genes, 
and directed links modeling gene interactions. 
Dynamical state of an individual gene can be de- 
scribed by a pair of variables representing con- 
centrations of protein and mRNA with station- 
ary time variations over the cell-cycle. The ar- 
chitecture of Gene Regulatory Networks con- 
tribute to cell's optimal functionality and dy- 
namical stability, which are both necessary for 
its survival. An interesting open question re- 
gards the potentials of the structure of a given 
Gene Regulatory Network to stabilize the dy- 
namics for a wider class of nonlinear dynami- 
cal systems, such as coupled chaotic maps. It 
has been recognized that chaotic maps on net- 
works of different types display a variety of col- 
lective effects including self-organized stable dy- 
namics and new types of attractors. In this pa- 
per we consider Gene Regulatory Network of 
bacterium Escherichia Coli and examine its abil- 
ity to induce coherent and stable dynamical pat- 
terns into the system of chaotic maps interacting 
with time delay through its directed regulatory 
links. 



1 Introduction 

The gene regulation is a process of fundamental impor- 
tance for the functioning and growth of biological cells. 
Gene regulatory system consists of genes which func- 
tion collectively by interacting through the transcrip- 
tion factors (activators or repressors), and produce the 
appropriate proteins in response to cell's needs pQ. 

A system of interacting genes can be seen as a net- 
work [2] , where nodes represent genes and directed links 
model the activatory and repressory interactions - such 
description is termed Gene Regulatory Network (GRN) 
[3]. Recent studies revealed various details of GRN's ar- 
chitecture, dynamics and functionality, and allowed ap- 
plications such as engineering of gene circuits [I] . GRN 
of many living organisms are known empirically, both 
in terms of chemistry of gene interactions and the de- 
tails of network topology j5j[6]. A well studied example 
is bacterium Escherichia Coli (E.Coli), where detailed 
experimental and theoretical investigations revealed all 
relevant gene interactions [BJ [7] . GRN typically involve 
different scales of functioning and have modular struc- 
ture, where groups of genes with specific connection 
patterns preform particular functions [U EJ [6J . 
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For modeling the interaction of biological units such 
as genes, one often uses simple dynamical systems that 
capture the essence of system's behavior over time [H 
Ell ED Q2] • In this way, the dynamics of gene regu- 
lations is modeled at various levels [T3] and using dif- 
ferent mathematical techniques [9]. The simplest ap- 
proaches include boolean networks [TU] and discrete- 
time maps [11], while more detailed analysis requires 
ID [H] or 2D [12] continuous-time ODE, that can also 
involve a time-delayed action [15]. Analytical studies 
of GRN models involving small networks revealed the 
complexity of their dynamical patterns [13], with the 
system of two interacting genes solved in detail [11} [T2] . 
General relationship between structure and function of 
large GRN was examined, with particular emphasis on 
their collective dynamics [16] . information processing 
[17j . flexibility |10j . and functional organization [7J. 

An important question concerning both small and 
large GRN structures is related to the stability of its 
emergent dynamics: GRN functioning ought to be very 
robust in order to assure the survival of the cell under 
various circumstances. How much stability and func- 
tionality can a certain network topology provide for a 
system of dynamical units coupled through its links? 
General relationship between topology and stability was 
studied for the case of static [18] and growing networks 
[19j . Stability of the GRN dynamics was examined from 
the Control theory prospective for noisy [20j and non- 
noisy case [15] . As time delay for all interactions is 
ubiquitous [21], stability and robustness of GRN were 
also investigated in the case of time- varying delays [22J , 
and asynchronous networks [23]. The dynamics at the 
"edge of chaos" characterized by zero Lyapunov ex- 
ponents was found in models of gene interaction [24j . 
Some recent studies [25] suggest importance of the crit- 
ical region between the order and chaos in the emergent 
dynamics of many biological networks. 

Networks of coupled maps are a useful paradigm 
for designing complex dynamical systems. Networks of 
maps have been extensively studied on different arti- 
ficial networks in the context of testing the effects in- 
duced by network topology on the emergent phenom- 
ena, where typically ID maps/ODE have been exam- 
ined [T71 [18]. Chaotic maps have been considered as 
suitable models that take into account certain features 
of gene interactions |llj . in addition to allowing the 
study of time delay effects on the emergent behaviors 
[21j . Nonlinear systems were used for modeling brain 
dynamics [22] and temporal fluctuations in gene regu- 
lations [23]. Coupled maps on empirical networks have 



been employed as simple tools for examining the ro- 
bustness of the network topology [26] . Recently, two- 
dimensional chaotic maps symmetrically coupled with 
time delay have been studied on artificial scale-free net- 
works. They are found to exhibit a rich variety of 
the emergent dynamical behaviors \27\ [28], including 
strange non-chaotic attractors [29] . 

The purpose of the present work is somewhat differ- 
ent: we focus on a realistic Gene Regulatory Network 
of bacterium E.Coli, and study its dynamical stability 
using a different class of nonlinear dynamical systems 
as our tool. In particular, we employ two-dimensional 
chaotic maps associated with the nodes, while vary- 
ing the coupling types and strengths. We determine 
the conditions under which the gene regulatory archi- 
tecture maintains its dynamical stability and describe 
the pathways towards the stable collective dynamical 
states, which contrasts the strongly chaotic dynamics 
exhibited by the isolated nodes. We describe the emer- 
gent states of the network by several quantitative mea- 
sures. 

For the nonlinear chaotic map associated with each 
node we use 2D Chirikov standard map [31], which 
is well understood and shares some formal properties 
(specifically the dimensionality of the phase space) with 
the actual gene dynamics, given at each node by the 
concentrations of protein and mRNA. It contains built- 
in periodicity in the angle variable, while being un- 
bounded in the action variable. We stress that this 
choice is primarily motivated by a tunable chaotic be- 
havior of the isolated standard map, which allows an 
easy quantification of the emergent non-chaotic dynam- 
ics that is a clear consequence of the network interac- 
tion among the maps. Our system is dissipative, despite 
involving conservative standard maps as units. Dissi- 
pativeness is a general feature of biological dynamics, 
which in our model arises as a consequence of the par- 
ticular coupling form that we chose, in addition to the 
network's directedness. 

The Network. We consider the version of E.Coli's 
GRN shown in Fig.[TJ which was empirically determined 
in 2003 [6] (data available in Rcf.[30j). Since we are in- 
terested in the induced collective dynamics, we consider 
only the largest connected component, which is com- 
posed of N = 328 nodes/genes (shown in Figs.E]&[8]). 
Smaller components have their own mutually unrelated 
dynamics. Compared to our previous works where the 
coupled maps were studied on non-directed networks 
[27] [28] , the focus here is on the effects that the direct- 
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Figure 1: Directed E.Coli's GRN with N = 423 genes 
marked by their biological names, according to the data 
from [61 [30] . Self-loops as well as the difference between 
activatory and repressory links are omitted. 



edness of links with time-delayed interactions have on 
the emergent collective dynamics and its stability. 

The paper is organized as follows: in Section [2j we 
introduce the system of chaotic maps coupled through 
the directed links of E.Coli's GRN and study briefly the 
case of SUM-rule coupling. Our main results are given 
in Section [3[ where we study the Spreading coupling 
between the maps and analyze in detail the collective 
dynamics in two characteristic regions of the coupling 
strengths. A short summary and the discussion of the 
results are given in Section [H 

2 Chaotic Maps on Directed E.Coli 
GRN 

In this Section we define the directed network of chaotic 
2D maps and introduce statistical and stability mea- 
sures that will be used to characterize its emergent dy- 
namics. We also show the main results for the case of 
SUM-rule coupling, which are compared to the results 
for non-directed case [27] . 

The Network of Chaotic Maps. We consider the 
largest connected component of the directed E.Coli's 
GRN, shown in Figs.[6]&[8l and assign to each node [i] = 
1, . . . N a Chirikov standard map which consists of two 



mutually coupled (action-angle) variables: 

/ x[i]' \ _ ( x[i] + y[i] + e sin(2-7rx[i]) modi \ 
V VW )~\ !/[*] +£sin(27nE[i]) J 

(1) 

representing a discrete version of the kicked rotor [31j . 
The chaotic parameter e is set to e = 0.9, thus assur- 
ing a strongly chaotic behavior of the isolated units. 
Each node is described by its dynamical state in two- 
dimensional phase space, represented by a point (x[i]t , y[i]t) 
at a discrete time t. The maps are coupled along the di- 
rected network links, with a coupling term that involves 
a one-step time delay in the coupled angle coordinates 
x between the network neighbors: 

( )=(l-a)( \ + 

(2) 

M(n, kf 1 ) defines the type of interaction between the 
nodes, which always depends on the coupling strength 
fi. Note that in a directed network /Q denotes only 
the neighbors having links towards the node [i] , and ac- 
cordingly, the interaction term M^fijkf 1 ) involves only 
the node's in-degree kf 1 . Some nodes may not have in- 
coming links {kf 1 = 0) and have thus vanishing inputs 
from other nodes. However, their motion is eventu- 
ally affected due to (1 — jj) term, which inhibits the 
chaotic diffusion |28j . Time delay is realized through 
the fact that the coupling term contains the node's up- 
dated value while the neighboring nodes' values 
(x[j]t) are not updated (for an update we here intend 
the action of the standard map denoted as ' and defined 
in Eq. {TJ). 

In this Section we briefly discuss the case with the 
coupling form 

M{^kT) = »k? (3) 

which resembles the SUM logic-gate GRN model: the 
input that a gene receives is the average of the inputs 
coming from all the neighbors [20] (a different type of 
coupling is discussed in Section [3l). It is important to 
stress that the neighboring nodes are coupled via their 
angular variables x only, while the dynamics of the ac- 
tion coordinate y is affected via standard map update 
Eq. ([1]), and by the prefactor (1 — /x) in Eq. ([2]). 

Types of Orbits and Statistical Measures Em- 
ployed. In the simulations with the dynamics defined 
by Eqs. ([T][2]), we start with a random selection of ini- 
tial conditions from the interval (x, y) € [0, 1] x [—1, 1] 
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for each GRN node, and a fixed coupling fx along the 
network links. The transients are taken as to = 10 5 
iterations. Two types of the emergent orbits are con- 
sidered: 

• orbit of the individual nodes: (x[i]t,y[i]t)t>t , [i] = 
1,...,N. 

• time-averaged orbit for a given node: (x[i], y[i\) = 
lim^oo Y?k=t ( x l i '\k,y[ i ]k) which reduces an 
entire orbit of a given node into a single phase 
space point. It qualitatively captures its motion 
after transients (or during the time-interval [to,t] 
if t < oo). 



exponents typically used with isolated chaotic systems 
[28] . FTMLE are a suitable measure for single nodes 
attached to the rest of the network, thus giving a nodc- 
by-node characterization of the network stability. 

• Return-times to Phase space Partitions represent an 
interesting measure of the non-ergodicity of the dynam- 
ics [27], and thus provides a suitable way of quantifying 
the collective effects. We monitor the time intervals 
between successive visits of an orbit to a phase space 
partition element. As discussed in Section [3J, for this 
purpose we consider 10 5 equal elements along the x- 
coordinate, while the cells remain unlimited along the 
y-coordinate. 



As mentioned above, we are testing the ability of the 
examined network structure to regularize the motion of 
the chaotic maps due to their mutual interactions along 
the directed links. The system's emergent dynamics will 
be characterized using following tools: 
• Periodic orbits for single nodes are defined to have a 
periodicity tt = t\ — to if for some t\ > to we have 



XI 



to 



x[i] tl 



XI 



to 



< 5 and 



vHto -l/[*]ti 



< 5 



Here we set 6 = 10~ 4 . Appearance of stable periodic 
orbits indicates a presence of regularity in the collective 
dynamics, in a clear contrast with the chaotic behavior 
of the isolated standard map. 

• Finite-time Maximal Lyapunov Exponent (FTMLE) 
for a point xo is defined as the maximal initial diver- 
gence rate between xo and the trajectories starting in 
its neighborhood J\f: 



A maa .(xo) = max ^ initial slope 



1 d(U T x, U t -k ) 
t d(x,xo) 



where the operator U T denotes the time-evolution ac- 
cording to Eqs. ([HE]) until the time r. The actual r- 
value used in this expression is determined for each par- 
ticular orbit, which can be identified as either strongly 
or weakly chaotic, or otherwise regular. We first deter- 
mine the character of the orbit by computing an approx- 
imate FTMLE using r = 10 3 iterations, and only then 
set the r-value used to compute the real FTMLE. Typ- 
ically, t ranges from 20 iterations for strongly chaotic 
to 200 iterations for stable orbits. For a given emergent 
orbit A^ ax -s are averaged for many initial points xo 
belonging to it, and )^ max = (A^ naa .(xo)) Xo is considered 
the characteristic FTMLE for this orbit, quantifying 
its stability. In opposition to the standard Lyapunov 



The Case of SUM-rule Coupling. We find that the 
GRN network with the SUM-rule coupling between the 
units given by Eq. ([2][3]) exhibits periodic orbits with dif- 
ferent periodicities for any non-zero //-value, as shown 
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Figure 2: 2D histogram of y- values in (a), and A' 
values in (b). The color-map shows the fraction of nodes 
having a certain value of y or ^ max as function of the 
coupling [i. Log-scale is used for the color-map, with the 
lowest (reference) value taken as 10~ 14 , corresponding 
to the deep blue color. 

in Fig. [2] In addition, self-organized non-periodic or- 
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bits are found at some nodes for the coupling strengths 
around u ~ 0.05. To compare with the case of non- 
directed (symmetrical) network from [27] , we study the 
clustering and stability of the emergent orbits. A his- 
togram of y[i] -values as function of fi is reported in 
Fig. [2K showing groups of nodes with clustered orbits. 
In Fig. [2b the histogram for the corresponding FTMLE 
values )t max is shown. Both histograms are obtained 
for a single set of initial conditions at every a- value and 
averaged over all the network nodes. 

These results are qualitatively similar to the case of 
non-directed networks: the system exhibits clusters of 
stable periodic orbits for any non-zero coupling, along 
with weakly chaotic orbits (small positive A^j 0a . ~ 0.1) 
around li ~ 0.05 for a some small fraction of nodes. 
Strongly chaotic orbits similar to the isolated standard 
map's ones are occurring for very small /i-values, and 
disappear for larger //-values. The number of clusters 
decreases and eventually shrinks to one. In analogy 
with the non-directed coupling, three dynamical regions 
can be distinguished: strongly chaotic (fi < 0.01), pe- 
riodic (0.01 < u < 0.04), and weakly chaotic with self- 
organized orbits (fj, > 0.04). The network dynamics for 
large coupling strengths (beyond the shown interval) re- 
mains stable with all orbits being periodic. We conclude 
that with the appropriate coupling form - SUM-rulc 
with the normalized inputs at each node - E.Coli's di- 
rected GRN is generally able to regularize the behavior 
of the coupled chaotic maps in a wide range of coupling 
strengths. 

3 Chaotic Maps on E.Coli GRN with 
Spreading Coupling 

We now modify our network of coupled maps by em- 
ploying the Spreading coupling form. In Eq. ([2]) we set: 



u = const. 



(4) 



We further test the potentials of the directed E.Coli's 
GRN to stabilize the behavior of chaotic standard map 
associated with its nodes, by investigating the emer- 
gent dynamics of Eqs. ([21&SD using the tools and the 
approach introduced in the previous Section. In con- 
trast to the SUM-rule coupling studied above, in this 
case the system is expected to be more sensitive to the 
chaotic behavior, given that better connected nodes re- 
ceive larger total inputs compared to the nodes with 
fewer links. With this type of coupling the dissipa- 
tiveness becomes unevenly distributed over nodes, with 



better linked nodes undergoing stronger dissipative ef- 
fects. As it will be shown below, this inhomogeneity 
of the dissipation dramatically affects the collective dy- 
namics. 

In Fig. [3] we report the fraction of non-periodic or- 
bits as function of li for all network nodes, averaged 
over many initial conditions. The profile of the curve 
at small /i-values resembles the one for the non-directed 
network and the SUM-rule cases: for all initial condi- 
tions, after a quick initial transient all nodes stabilize 
into different periodic orbits. However, with the in- 
crease of the coupling strength a fraction of nodes is 
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Figure 3: Fraction of non-periodic orbits as function of 
the coupling strength fx, for all nodes and averaged over 
many initial conditions. Inset: zoom to the dynamical 
transition region around \x ~ 0.022. 



found around u ~ 0.02, whose emergent orbits remain 
non-periodic for all initial conditions (inset in Fig. [3]). 
Generally, with further increase of /x the fraction of 
nodes with non-periodic orbits fluctuates, but overall 
slowly increases. Properties of these nodes and their 
orbits will be studied in detail below. 

The 2D histogram of -values, shown in Fig. [4^, 
is reflecting different dynamical regions with the oc- 
currence of the non-periodic orbits, in agreement with 
Fig. El For small /x-values we find a fully clustered or- 
ganization of nodes with periodic orbits only, whereas 
for u > 0.02 only a few separate clusters of orbits can 
be identified. Contrary to the SUM-rule case, here we 
find that for fx > 0.08 the network displays mainly non- 
clustered orbits, co-existing with a few remaining clus- 
ters. The situation is also reflected in the 2D histogram 
of the FTMLE A^^, shown in Fig. [4)3. For small non- 
zero couplings li, non-positive A^ na3 , can be found at 
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Figure 4: 2D histogram of y[i]-values in (a) and A 
values in (b) for directed network dynamics Eqs. ([2]&;[3D, 
computed for a single initial condition and for all nodes. 
The log-scale color-map is as explained in Fig.[2j 



all nodes. However, destabilization of nodes occurs for 
fx > 0.02, where we find a spectrum of positive X max - 
values. In this region of coupling strengths X max > 1 is 
found at some nodes, indicating their strongly chaotic 
behavior, similarly to the isolated standard map. Apart 
from exhibiting the unstable dynamics with positive 
FTMLE at the majority of nodes, the 2D histogram 
of X max -values for fj, > 0.022 exhibits a characteristic 
pattern, as visible in Fig.Hb, that includes a spectrum 
of small FTMLE X max > 0. In addition, this patterns 

> 

max r*j " 



suggests that the weakly chaotic nodes having A* 
occur as a large population of the network nodes, in 
opposition to a much smaller fraction of nodes with 
X max > 1. In what follows we study in detail the emer- 
gent dynamics at all nodes on the network for two par- 
ticular coupling strengths: at the onset of destabiliza- 
tion for [i ~ 0.022, and within the unstable region at 
\i = 0.05. 

Dynamics at the Onset of Destabilization \i ~ 
0.022. The collective dynamics of the directed E.Coli's 



GRN described by Eq. ([2]) with the Spreading coupling 
rule Eq. (j4]) becomes partially unstable for [i > 0.02, 
exhibiting a mixture of stable and weakly chaotic or- 
bits. As shown in the inset in Fig.[3j the destabilization 
occurs gradually, resembling a continuous phase tran- 
sition. For the coupling strengths close to the thresh- 
old [i ~ 0.022 we examine the fraction of non-periodic 
orbits at each node for many initial conditions. The 
results are shown in Fig. [5^. We can identify the nodes 
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Figure 5: For each network node [i] the color-map shows 
the fraction of initial conditions that: (a) lead to non- 
periodic orbits as a function of /i in the vicinity of the 
threshold, and (b) lead to a given A^ aa ,-exponent for 
the threshold coupling fj, = 0.022. 



that gradually lose periodic orbits as the intensity of 
coupling increases beyond fi = 0.022. These nodes seem 
to undergo the process almost simultaneously, while the 
rest of the network appears to be unaffected. In order 
to further analyze the nature of the non-periodic orbits, 
we compute and show in Fig. [5b the histogram of the 
A^j a3 .-values for each network node separately, at this 
particular coupling strength of \i = 0.022. Indeed, the 
nodes that lose periodic orbits display A^ aa , > 0. Posi- 
tive A^ aa . are smaller than 1, indicating weakly chaotic 
effects at the unstable nodes. It is striking feature of 
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this directed network that the instability develops on a 
specific group of nodes, rather than over the whole net- 
work. Moreover, this instability remains confined to the 
initially destabilized sub-network even after the transi- 
tion at \x = 0.022. Whereas, the rest of the network 
manages to maintain its regular stable dynamics. 

In Fig. [6] we show the largest connected component 
of the directed E.Coli's GRN indicating the node stabil- 
ity: nodes with periodic orbits at \i = 0.022 are shown 
in blue and non-periodic ones in red color. The unstable 
nodes are mutually connected and form a sub-network 
composed of 24 nodes, and include the hub-node (a 
list of biological names of these unstable nodes / genes is 
given in Ref . |32] ) . The unstable sub-network has both 




Figure 6: Largest connected component of directed 
E.Coli GRN. Nodes exhibiting non-periodic orbits at 
H = 0.022 are marked by red, and those with periodic 
orbits by blue color. 

incoming and outgoing links with the rest of the net- 
work. This suggests a non-trivial interplay between the 
collective dynamics and the network architecture, re- 
sponsible for keeping the unstable dynamics co-existing 
with the stable one within the same network. We run 
the simulations for longer than 10 6 iterations to confirm 
the stationarity of these results. 

In Fig. [T^i, b &: c three typical attr actors/orbits ap- 
pearing on single unstable nodes at [i = 0.022 are shown, 
exhibiting different structural patterns within reduced 
of the phase space. The hub-node always displays 
the strange attractor from Fig. [7^ with a fractal struc- 
ture, regardless of the initial conditions (A^ aa . ~ 0.4, 
fractal dimension df ~ 1.4). We examine the statisti- 
cal properties of these single-node orbits by computing 



their return time distributions with respect to the phase 
space partitions, defined above. The obtained distribu- 
tions are fitted with ^-exponential function defined as 
[28]: 

e q (x) = B q (l-(l- q )^y~ q (5) 

The results for all three orbits are shown in Fig.[7]i. 
The distributions indicate long-time correlations in the 
return times, compatible with the power-law tail hav- 
ing q > 1. The trajectories explore different parts of 
the phase space non-uniformly, in a clear opposition to 
the exponential distribution (for q = 1), which charac- 
terizes full ergodicity of the chaotic dynamics exhibited 
by the isolated nodes. Occurrence of the power-law 
tails of the distributions suggests the presence of self- 
organization effects, arising as a consequence of the in- 
teraction among the nodes, which leads to the emergent 
non-periodic orbits with A > 0. 

In a sharp contrast with what was observed in the 
previous Section, the directed network with the dynam- 
ics given by Eqs. (|2]&;|3]) remains destabilized above the 
coupling intensity of /x = 0.022. Also, contrary to the 
uniformity of the emergent states in the case of the sym- 
metrically coupled maps, in the present case we find 
that several patterns of different dynamical nature can 
be simultaneously present in different parts of the net- 
work. 

Collective Dynamics inside Chaotic Region at 
/i = 0.05. In the region of coupling strength fi > 0.022 
the GRN nodes are exhibiting three types of behav- 
ior: stable periodic orbits, weakly chaotic orbits, and 
strongly chaotic orbits (cf. Figs.[3]&[4]). In Fig. [8] we 
show the state of the network at fj, = 0.05, where the 
non-periodic nodes are indicated by red color. The un- 
stable behavior is again localized to a connected sub- 
network, which is integrated within the rest of the net- 
work. The unstable sub-network comprises 10% of the 
total number of nodes and includes the unstable sub- 
network found in case of fi = 0.022 (with additional 
10 nodes listed in Ref. [33]). Despite stronger coupling, 
the network still manages to contain its dynamically 
unstable part within a small sub-network, allowing for 
regular stable behavior of the rest of the nodes. 

In Fig. [9^, we show 2D histogram of y [i]-values nodc- 
by-node at \x = 0.05. While most of the nodes' orbits 
are clustered into three main clusters for all initial con- 
ditions (horizontal lines in this plot), some other nodes 
have a wide range of possible y[i]-values away from a 
clear clustering pattern. For comparison, in Fig. [9}} 
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Figure 7: Examples of typical orbits for unstable single nodes at [i = 0.022: (a) the attractor of the hub-node 
2 (gene biological name "acrAB"), (b) an orbit of the node 214 (gene "kbl-tdh"), and (c) an orbit of the node 
46 (gene "aspV"). Only representative parts of the orbits are shown, (d) distributions of return times to phase 
space partition for these single-node orbits, fitted with the g-exponential expression Eq. ([5]), with the respective 
q- values listed. Phase space is partitioned into 10 5 equally spaced cells in x-coordinate (y-coordinate unbounded). 



we also show the histogram of A^ aa .-values for each 
node on the network and the fixed coupling strength 
[i = 0.05. As the Fig. [9b shows, the majority of nodes 
appear to have stable dynamics (A^ aa . < 0). For the 
unstable nodes we find the values of )J max to vary from 
Knax ~ 0.1 (also present at the instability threshold 
fx = 0.022), to a strong chaos with A^ ax > 1. We also 
find that some of these nodes alternate between stable 
and weakly unstable dynamics. In our coupled maps 
dynamics on directed GRN, each node attains a certain 
spectrum of "roles" in the network emergent behavior, 
depending on the particular choice of the initial con- 
ditions and the connection patterns, thus implying the 
roles of other nodes connected to it. 

Contrary to the case of /i = 0.022 shown in Fig.0 at 
the strong coupling the hub node does not display any 
specific attractor, but it shows strongly chaotic motion 
regardless of the initial conditions. However, the behav- 



ior of some other non-periodic nodes shows a variety of 
patterns that are structurally and statistically robust 
to the initial conditions. We illustrate this in Fig.llO[ 
where three typical emergent orbits of two selected un- 
stable nodes are shown, along with the corresponding 
distributions of the return times. The key geometrical 
and statistical properties of the orbits seem to persist 
for varying initial conditions and are related with the 
observed fluctuations in y[i]-values. The shapes of re- 
turn times distributions for these orbits are shown in 
Fig.llOb &: d. Again, these distributions can be fitted 
with the expression Eq. ([5]) with q > 1, indicating an 
self-organized collective dynamics, away from the fully 
chaotic dynamics of the isolated maps. 

The emergent dynamics of our coupled chaotic maps 
system possesses a certain flexibility with respect to the 
initial conditions (which can be seen as the environ- 
mental inputs to the cellular GRN). It is marked by 
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Figure 8: Coupled chaotic maps on the directed E.Coli 
GRN for \i = 0.05: non-periodic nodes are localized 
in the red sub-network, while the rest of nodes (blue) 
manage to attain the periodic orbits. 



the finite range of possible y[i]-values and A^ ax -values. 
On the other hand, the statistical properties of the dy- 
namics, i.e., the return-time distributions of the angular 
variable, for particular unstable nodes/genes, appear to 
be quite robust to variations of the initial conditions. 
Features of this type might be related to the biological 
origin of the studied network: the network's architec- 
ture provides certain level of functional adaptability and 
robustness for operation of its units - genes, even when 
the inherently chaotic dynamics at its units is present. 

4 Discussion and Conclusions 

We have studied the emergent dynamics and stability 
of two-dimensional chaotic standard maps coupled with 
time delay along the directed regulatory links of the 
largest connected component of Escherichia Coli's Gene 
Regulatory Network. The coupling with SUM-rule and 
Spreading-rule between the connected units have been 
considered, and the emergent dynamical properties were 
investigated through non-periodic orbits, and by using 
Finite-time Maximal Lyapunov Exponents y^ max and 
distributions of return-times to the phase space par- 
titions. We have demonstrated that the examined em- 
pirical network structure is able to induce a coherent 
collective dynamics of the coupled chaotic maps. The 
emergent behavior appears to depend primarily on the 
network structure, but also on the type and the strength 



^ 
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node [i] 

(b) 

Figure 9: For each node [i] the color-map shows the 
fraction of initial conditions for \i = 0.05 leading to: a 
given value of A^ aa ,-exponent in (a) and a given y[i\- 
value in (b). Log-scale is used as with Fig.[2j 



of coupling. 

Our numerical results suggest that different mech- 
anisms are in a combined manner contributing to the 
observed collective dynamical behavior of the coupled 
maps on network: 

• Self- organization among dynamically coupled units 
occurs over time, where the state of each unit ad- 
justs to the states of its neighbors in various ways, 
depending on the couplings and the local dynam- 
ics of its two variables, without any external in- 
teraction involved. 

• Anomalous diffusion in our coupled map network 
is demonstrated for the trajectories of some nodes, 
which are shown to selectively fill the phase space. 

• Non-ergodicity, related with the anomalous diffu- 
sion and the self-organized behavior of the cou- 
pled maps, is a general property of complex dy- 
namical systems. In the present case it might be 
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Figure 10: Examples of three emergent orbits (marked by different colors) for single nodes at \i = 0.05: the node 
86 (gene "zraP") in (a), and the node 119 (gene "fdhF") in (b). The distributions of return times to phase space 
partitions are shown for the orbits from (a) in (c), and for the orbits from (b) in (d). Phase space partitioning 
is done identically to the case from the Fig. UJ 



caused by the underlying network structure, com- 
bined with the dissipatively coupled maps on it. 

• Synchronization occurring between different nodes 
is documented by the appearance of the clusters of 
nodes with a well defined y [f]-values (cf. Figs.[3]&;[9j) 
Mutually synchronized nodes tend to form a char- 
acteristic pattern over the network which we have 
not studied for this case (results for non-directed 
networks can be found in J27J [28] ) . 

• Phenomena at the edge- of- chaos or a dynamical 
phase-transition are often reported (see e.g. |24[ 
[23] ) in the extended dynamical systems when the 
Lyapunov exponents vanish. Our results also sug- 
gest that the attractors with vanishingly small 
finite-time Lyapunov exponents A^ naa; might play 
an important role in the occurrence and the struc- 
ture of the collective dynamical states on the stud- 
ied GRN. 



The structural properties of the network are closely 
related with the observed dynamical features and their 
stability. The flexibility of the emergent orbits to the 
initial conditions indicates the network's ability to effi- 
ciently respond to the environmental inputs by adapt- 
ing dynamics of its units accordingly. Also, the results 
for Spreading coupling suggest the existence of an op- 
timal range of interaction strength where the network 
maintains its dynamical stability. The fact that the 
instability appearing at strong inter-node coupling re- 
mains localized to a small sub-network indicates the 
set of nodes/links which are crucial for the stability of 
this network structure. These nodes/links may be also 
used as a target for dynamical manipulation of the net- 
work. The unstable sub-network is dynamically linked 
to the stable part of the network, suggesting a complex 
balance of the global network's behavior, that needs 
additional study. 

We hope that our work traces the appropriate method- 



ic) 



ology for the study of collective dynamics in coupled 
multi-dimensional chaotic maps with time delay on em- 
pirical networks, in particular in the context of search- 
ing for a network architecture with robust and dynam- 
ically controllable behavior. 
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